function generateMfuncResidualEquation_file(allModelParams,x,y,xp,yp,nx1,ne,h11,heteroShock,strucOfArrays,positionFile,nameOfFunction)
% This function generates an M function which is the objective function
% when solving the model by projection
% The results are saved in the file ResidualEquationProj.m

% We start deleting the old version of the file - if it exists
%nameOfFunction = 'ResidualEquationProj.m';
saveFile = [positionFile,'\',nameOfFunction];
if exist(saveFile,'file') > 0
    delete(saveFile)
end

text = ['function [res,resEuler] = ',nameOfFunction(1:end-2),'(thetaVec,setup)'];
dlmwrite(saveFile,text,'-append','delimiter','');
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '%% Getting some information from setup ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'scaleX_in_g = setup.scaleX_in_g;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'x_cu        = setup.xGrid;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'numX        = size(x_cu,1);';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'nNodes      = setup.n_nodes;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'wNodes      = transpose(setup.weight_nodes);';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'epsNodes    = setup.epsi_nodes;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'params      = setup.params;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'appOrder    = setup.appOrder;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'nx          = setup.nx;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'ny          = setup.ny;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'nx1         = setup.nx1;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'gSS         = setup.gSS;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'hSS         = setup.hSS;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '% Unfolding thetaVec';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'thetaProj   = zeros(setup.nDim1Theta,setup.nDim2Theta);';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'countIdx = 0;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'for i=1:setup.nDim1Theta';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '    nn = countIdx + sum(setup.select(i,:));';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '    thetaProj(i,setup.select(i,:) == 1) = thetaVec(countIdx+1:nn);';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '    countIdx = nn;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'end';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'if setup.endoMeanGrid == 1';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '    % Re-centering the xGrid using the projection solution';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '    meanX = meanProj(thetaProj,setup);';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '    x_cu = x_cu+repmat(transpose(meanX),size(x_cu,1),1);';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'end';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'g0          = thetaProj(:,1); ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'gx          = thetaProj(:,2:1+nx); ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'if appOrder > 1 ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '    gxxV    = thetaProj(:,nx+2:1+nx+nx*(nx+1)/2); ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'end ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'if appOrder > 2 ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '    gxxxV   = thetaProj(:,1+nx+nx*(nx+1)/2+1:1+nx+nx*(nx+1)/2+nx*(nx+1)*(nx+2)/6); ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'end ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '%Unfold params ';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:length(allModelParams)
    text = [allModelParams{i}, '= params.',allModelParams{i},';'];
    dlmwrite(saveFile,text,'-append','delimiter','');
end 
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '%% Building the residuals ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '% Getting the states at time t ';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:length(x)
     text = [char(x(i)),' = hSS(',num2str(i),',1) + repmat(x_cu(:,', num2str(i),'),1,nNodes);'];
     dlmwrite(saveFile,text,'-append','delimiter','');
end
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '% Getting the controls at time t and part of the states at time t+1  ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '% yVec_cu: numX * ny ';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:length(x)
    text = ['x',num2str(i),'_cu ',' = x_cu(:,',num2str(i),')/scaleX_in_g(',num2str(i),',1);'];
    dlmwrite(saveFile,text,'-append','delimiter','');
end
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:length(y)
    name = [char(y(1,i)),' = ','g0(',num2str(i),',1)'];
    for j=1:length(x)
        name = [name,['+gx(',num2str(i),',',num2str(j),').*x',num2str(j),'_cu']];
    end
    text = [name,';'];
    dlmwrite(saveFile,text,'-append','delimiter','');
end
text = 'if appOrder > 1';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:length(y)
    name = [char(y(1,i)),' = ',char(y(1,i))];
    idx = 0;
    for alfa1=1:length(x)
        for alfa2=alfa1:length(x)
            idx = idx + 1;
            if alfa1 == alfa2
                name = [name,['+gxxV(',num2str(i),',',num2str(idx),').*x',num2str(alfa1),'_cu','.*x',num2str(alfa2),'_cu']];
            else
                name = [name,['+2.*gxxV(',num2str(i),',',num2str(idx),').*x',num2str(alfa1),'_cu','.*x',num2str(alfa2),'_cu']];
            end
        end
    end
    text = ['    ',name,';'];
    dlmwrite(saveFile,text,'-append','delimiter','');
end
text = 'end';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'if appOrder > 2';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:length(y)
    name = [char(y(1,i)),' = ',char(y(1,i))];
    idx = 0;
    for alfa1=1:length(x)
        for alfa2=alfa1:length(x)
            for alfa3=alfa2:length(x)
                idx = idx + 1;
                if alfa1 == alfa2 && alfa2 == alfa3
                    name = [name,['+gxxxV(',num2str(i),',',num2str(idx),').*x',num2str(alfa1),'_cu','.*x',num2str(alfa2),'_cu','.*x',num2str(alfa3),'_cu']];
                elseif alfa1 == alfa2 && alfa2 ~= alfa3
                    name = [name,['+3.*gxxxV(',num2str(i),',',num2str(idx),').*x',num2str(alfa1),'_cu','.*x',num2str(alfa2),'_cu','.*x',num2str(alfa3),'_cu']];
                elseif alfa1 ~= alfa2 && alfa2 == alfa3
                    name = [name,['+3.*gxxxV(',num2str(i),',',num2str(idx),').*x',num2str(alfa1),'_cu','.*x',num2str(alfa2),'_cu','.*x',num2str(alfa3),'_cu']];
                elseif alfa1 ~= alfa2 && alfa2 ~= alfa3
                    name = [name,['+6.*gxxxV(',num2str(i),',',num2str(idx),').*x',num2str(alfa1),'_cu','.*x',num2str(alfa2),'_cu','.*x',num2str(alfa3),'_cu']];
                else
                    [alfa1,alfa2,alfa3];
                end
            end
        end
    end
    text = ['    ',name,';'];
    dlmwrite(saveFile,text,'-append','delimiter','');    
end
text = 'end';
dlmwrite(saveFile,text,'-append','delimiter','');
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'xVec_cup = zeros(numX,nx);';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'mx = setup.mx;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'if mx > 0 ';
dlmwrite(saveFile,text,'-append','delimiter','');
nx11 = length(h11);
for i=1:nx11
    TEMP = char(h11(i,:));
    TEMP = strrep(TEMP,'*','.*');
    TEMP = strrep(TEMP,'/','./');
    TEMP = strrep(TEMP,'^','.^');
    text = ['    ',char(xp(1,i)),' = ',TEMP,';'];
    dlmwrite(saveFile,text,'-append','delimiter','');
end
text = 'end ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'if setup.myx > 0 ';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:nx1-nx11
    text = ['    xVec_cup(:,mx+',num2str(i),') = ',char(y(1,i)),'(:,1)-gSS(',num2str(i),',1);'];
    dlmwrite(saveFile,text,'-append','delimiter','');
end
text = 'end ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'diaghx = diag(setup.hx);';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'xVec_cup(:,nx1+1:nx)  = x_cu(:,nx1+1:nx).*repmat(transpose(diaghx(nx1+1:nx)),numX,1); ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '% The controls ';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:length(y)
    name = [char(y(1,i)),' = repmat(',char(y(1,i)),',1,nNodes);'];
    text = name;
    dlmwrite(saveFile,text,'-append','delimiter','');
end
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '%Adding shocks to the states at t + 1 ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'diagEta = diag(setup.eta(nx1+1:nx,:)); ';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:nx1
    name = ['x',num2str(i),'_cup  = repmat(xVec_cup(:,',num2str(i),'),1,nNodes);'];
    text = name;
    dlmwrite(saveFile,text,'-append','delimiter','');
end
for i=1:ne
    name = ['x',num2str(nx1+i),'_cup  = repmat(xVec_cup(:,',num2str(nx1+i),'),1,nNodes) +',...
         '2./(exp(-',char(heteroShock(i,1)),'.*x_cu(:,',num2str(nx1+i),'))+1).*diagEta(',...
         num2str(i),',1).*repmat(epsNodes(',num2str(i),',:),numX,1);'];
    text = name;
    dlmwrite(saveFile,text,'-append','delimiter','');
end
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '% Scaling the states at time t+1';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:length(x)
    name = ['x',num2str(i),'_cup  = x',num2str(i),'_cup/scaleX_in_g(',num2str(i),',1);'];
    text = name;
    dlmwrite(saveFile,text,'-append','delimiter','');
end
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '% Getting the controls at time t+1; each element is numX * nNodes';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:length(yp)
    name = [char(yp(1,i)),' = ','g0(',num2str(i),',1)'];
    for j=1:length(x)
        name = [name,['+gx(',num2str(i),',',num2str(j),').*x',num2str(j),'_cup']];
    end
    text = [name,';'];
    dlmwrite(saveFile,text,'-append','delimiter','');
end
text = 'if appOrder > 1';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:length(yp)
    name = [char(yp(1,i)),' = ',char(yp(1,i))];
    idx = 0;
    for alfa1=1:length(x)
        for alfa2=alfa1:length(x)
            idx = idx + 1;
            if alfa1 == alfa2
                name = [name,['+gxxV(',num2str(i),',',num2str(idx),').*x',num2str(alfa1),'_cup','.*x',num2str(alfa2),'_cup']];
            else
                name = [name,['+2.*gxxV(',num2str(i),',',num2str(idx),').*x',num2str(alfa1),'_cup','.*x',num2str(alfa2),'_cup']];
            end
        end
    end
    text = ['    ',name,';'];
    dlmwrite(saveFile,text,'-append','delimiter','');
end
text = 'end';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'if appOrder > 2';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:length(y)
    name = [char(yp(1,i)),' = ',char(yp(1,i))];
    idx = 0;
    for alfa1=1:length(x)
        for alfa2=alfa1:length(x)
            for alfa3=alfa2:length(x)
                idx = idx + 1;
                if alfa1 == alfa2 && alfa2 == alfa3
                    name = [name,['+gxxxV(',num2str(i),',',num2str(idx),').*x',num2str(alfa1),'_cup','.*x',num2str(alfa2),'_cup','.*x',num2str(alfa3),'_cup']];
                elseif alfa1 == alfa2 && alfa2 ~= alfa3
                    name = [name,['+3.*gxxxV(',num2str(i),',',num2str(idx),').*x',num2str(alfa1),'_cup','.*x',num2str(alfa2),'_cup','.*x',num2str(alfa3),'_cup']];
                elseif alfa1 ~= alfa2 && alfa2 == alfa3
                    name = [name,['+3.*gxxxV(',num2str(i),',',num2str(idx),').*x',num2str(alfa1),'_cup','.*x',num2str(alfa2),'_cup','.*x',num2str(alfa3),'_cup']];
                elseif alfa1 ~= alfa2 && alfa2 ~= alfa3
                    name = [name,['+6.*gxxxV(',num2str(i),',',num2str(idx),').*x',num2str(alfa1),'_cup','.*x',num2str(alfa2),'_cup','.*x',num2str(alfa3),'_cup']];
                else
                    [alfa1,alfa2,alfa3];
                end
            end
        end
    end
    text = ['    ',name,';'];
    dlmwrite(saveFile,text,'-append','delimiter','');    
end
text = 'end';
dlmwrite(saveFile,text,'-append','delimiter','');
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '% The states at time t + 1 ';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:length(x)
     text = [char(xp(i)),' = hSS(',num2str(i),',1) + x',num2str(i),'_cup*scaleX_in_g(',num2str(i),',1);'];
     dlmwrite(saveFile,text,'-append','delimiter','');
end
text = '% The Euler errors ';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:length(y)
    TEMP = char(strucOfArrays.('f')(i,:));
    TEMP = strrep(TEMP,'*','.*');
    TEMP = strrep(TEMP,'/','./');
    TEMP = strrep(TEMP,'^','.^');
    text = ['f',num2str(i),' = ', TEMP,';'];
    dlmwrite(saveFile,text,'-append','delimiter','');    
end
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:length(y)
    text = ['f',num2str(i),' = f',num2str(i),'*wNodes;'];
    dlmwrite(saveFile,text,'-append','delimiter','');
end
text = '% Numerical integration ';
dlmwrite(saveFile,text,'-append','delimiter','');
name = ['['];
for i=1:length(y)
    if i < length(y)
        name = [name,['f',num2str(i),',']];
    else
        name = [name,['f',num2str(i)]];
    end
end
text = ['resEuler = ',name,'];'];
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'resEulerWeight = resEuler.*(1./(1+sum(x_cu.^2,2)));';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'res = resEulerWeight(:);';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'end';
dlmwrite(saveFile,text,'-append','delimiter','');
end

